home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / LIBRARY / BPL70N16 / ARISOURC.ZIP / FPSRT.ASM < prev    next >
Assembly Source File  |  1993-03-07  |  10KB  |  198 lines
  1.  
  2. ; *******************************************************
  3. ; *                                                     *
  4. ; *     Turbo Pascal Runtime Library Version 7.0        *
  5. ; *     Real Square Root                                *
  6. ; *                                                     *
  7. ; *     Copyright (C) 1989-1993 Norbert Juffa           *
  8. ; *                                                     *
  9. ; *******************************************************
  10.  
  11.              TITLE   FPSRT
  12.  
  13.  
  14. CODE         SEGMENT BYTE PUBLIC
  15.  
  16.              ASSUME  CS:CODE
  17.  
  18. ; Externals
  19.  
  20.              EXTRN   HaltError:NEAR
  21.  
  22. ; Publics
  23.  
  24.              PUBLIC  RSqrt
  25.  
  26. ;-------------------------------------------------------------------------------
  27. ; RSqrt computes the square root of it argument. For a negative argument,
  28. ; runtime error 207 is invoked. The routine uses a estimate-and-correct method
  29. ; similar to that used in the division routine. This algorithm is based on the
  30. ; longhand method taught in school. Since there can be no tie case in rounding
  31. ; when computing the square root, no guard and sticky flags are needed to round
  32. ; to nearest or even in compliance with the IEEE floating point standard.
  33. ;
  34. ; INPUT:     DX:BX:AX  argument
  35. ;
  36. ; OUTPUT:    DX:BX:AX  square root of argument
  37. ;
  38. ; DESTROYS:  AX,BX,CX,DX,SI,DI,Flags
  39. ;-------------------------------------------------------------------------------
  40.  
  41. ROOT_EXT     PROC    FAR
  42. $zero_sqrt:  RET                       ; return argument
  43. $sqrt_err:   MOV     AX, 0CFh          ; load error code
  44.              JMP     HaltError         ; execute error handler
  45. ROOT_EXT     ENDP
  46.  
  47.              ALIGN   4
  48.  
  49. RSqrt        PROC    FAR
  50.              OR      AL, AL            ; argument zero ?
  51.              JZ      $zero_sqrt        ; yes, return zero
  52.              OR      DH, DH            ; argument negative ?
  53.              JS      $sqrt_err         ; yes, error
  54.              PUSH    BP                ; save TURBO-Pascal frame pointer
  55.              SHR     AL, 1             ; divide biased exponent by 2
  56.              SBB     CL, CL            ; CL = 0, if expo even, else CL = -1
  57.              ADC     AL, 40h           ; make new biased exponent
  58.              PUSH    AX                ; save new exponent
  59.              XOR     AL, AL            ; clear LSB of a3
  60.              OR      DH, 80h           ; set hidden bit
  61.              NEG     CL                ; CL = 0, if expo even, else CL = 1
  62.              SHR     DX, CL            ; divide argument
  63.              RCR     BX, CL            ;  by 2 if
  64.              RCR     AX, CL            ;   expo odd
  65.              XCHG    AX, CX            ; argument in DX:BX:CX
  66.              MOV     SI, DX            ; save a1
  67.              MOV     DI, BX            ; save a2
  68.              MOV     BP, CX            ; save a3
  69.              MOV     BH, DH            ; get MSB of a1
  70.              STC                       ; generate 4 bit approximation
  71.              RCR     BH, 1             ;  in hi-nibble of BH
  72.              NEG     AL                ; adjust approximation (subtract 10)
  73.              AND     AL, 10            ;  for arguments
  74.              SUB     BH, AL            ;   between 4000 0000 00 and 7FFF FFFF FF
  75.              MOV     AL, 0FFh          ; default quotient = FFh
  76.              CMP     DH, BH            ; division overflow ?
  77.              JNB     $no_div0          ; yes, use default quotient
  78.              MOV     AX, DX            ; get a1
  79.              DIV     BH                ; compute a1 / approximation
  80. $no_div0:    ADD     BH, AL            ; add result to approximation
  81.              RCR     BH, 1             ;  and divide by 2 yields 8 bit approx.
  82.              MOV     BL, 0FFh          ; BX >= Trunc(Sqrt(a1))
  83.              MOV     AX, 0FFFFh        ; default quotient
  84.              CMP     DX, BX            ; division overflow ?
  85.              JNB     $no_div1          ; yes, use default quotient
  86.              MOV     AX, DI            ; get a2
  87.              DIV     BX                ; compute a1:a2 / approximation
  88. $no_div1:    ADD     AX, BX            ; add result to approximation
  89.              RCR     AX, 1             ;  and divide by 2 yields 16 bit approx.
  90.              MOV     BX, AX            ; save q1
  91.              MUL     AX                ; DX:AX = q1*q1
  92.              SUB     DI, AX            ; compute
  93.              SBB     SI, DX            ;  remainder (SI:DI)
  94.              JNC     $quot1_ok         ; remainder must be positive
  95.  
  96.              ALIGN   4
  97.  
  98. $corr_quot1: DEC     BX                ; try next smaller quotient q1
  99.              STC                       ; adjust
  100.              ADC     DI, BX            ;  remainder
  101.              ADC     SI, 0             ;   to new
  102.              ADD     DI, BX            ;     value
  103.              ADC     SI, 0             ;      of quotient
  104.              JS      $corr_quot1       ; until remainder becomes positive
  105. $quot1_ok:   XCHG    AX, CX            ; concat rest of argument to remainder
  106.              MOV     DX, DI            ; get remainder lo-word
  107.              SHR     SI, 1             ; divide remainder (from here SI=0 !!!)
  108.              RCR     DX, 1             ;  by two so it fits into 32 bits
  109.              RCR     AX, 1             ; DX:AX = remainder / 2
  110.              MOV     CX, 0FFFFh        ; load default quotient
  111.              CMP     DX, BX            ; division overflow ?
  112.              JNB     $no_div2          ; yes, skip division and use default quot
  113.              DIV     BX                ; AX = q2
  114.              XCHG    AX, CX            ; BX:CX = q1:q2
  115. $no_div2:    MOV     AX, CX            ; get q2
  116.              MUL     BX                ; q1*q2
  117.              SUB     BP, AX            ; subtract product
  118.              SBB     DI, DX            ;  from
  119.              SUB     BP, AX            ;   old
  120.              SBB     DI, DX            ;    remainder
  121.              MOV     AX, CX            ; get q2
  122.              MUL     AX                ; q2*q2
  123.              NEG     AX                ; compute
  124.              SBB     BP, DX            ;   new
  125.              SBB     DI, SI            ;    remainder (DI:BP:AX)
  126.              JNS     $quot2_ok         ; remainder must be positive
  127. $corr_quot2: DEC     CX                ; try next smaller value for q2
  128.              STC                       ; adjust
  129.              ADC     AX, CX            ;  value
  130.              ADC     BP, BX            ;   of remainder
  131.              ADC     DI, SI            ;    according
  132.              ADD     AX, CX            ;     to changed
  133.              ADC     BP, BX            ;      value
  134.              ADC     DI, SI            ;       of q2
  135.              JS      $corr_quot2       ; until remainder positive
  136. $quot2_ok:   MOV     SI, AX            ; DI:BP:SI = remainder, BX:CX = quot
  137.              MOV     DX, BP            ; divide
  138.              SHR     DI, 1             ;  remainder by two and
  139.              RCR     DX, 1             ;   stuff 32 most significant bits
  140.              RCR     AX, 1             ;    into DX:AX
  141.              MOV     DI, 0FFFFh        ; load default quotient q3
  142.              CMP     DX, BX            ; division overflow ?
  143.              JNB     $no_div3          ; yes, use default quotient and skip div.
  144.              DIV     BX                ; AX = q3
  145.              XCHG    AX, DI            ; BX:CX:DI = q1:q2:q3, BP:SI:0:0 = rem
  146. $no_div3:    MOV     AX, DI            ; get q3
  147.              MUL     BX                ; q1*q3
  148.              SUB     SI, AX            ; subtract
  149.              SBB     BP, DX            ;  2*q1*q3
  150.              SUB     SI, AX            ;   from
  151.              SBB     BP, DX            ;    remainder
  152.              MOV     AX, DI            ; get q3
  153.              MUL     CX                ; q2*q3
  154.              PUSH    BX                ; save q1
  155.              XOR     BX, BX            ; BP:SI:BX:0 = remainder
  156.              SUB     BX, AX            ; subtract
  157.              SBB     SI, DX            ;  2*q2*q3
  158.              SBB     BP, 0             ;   from
  159.              SUB     BX, AX            ;    remainder
  160.              SBB     SI, DX            ;     "
  161.              SBB     BP, 0             ;      "
  162.              MOV     AX, DI            ; get q3
  163.              MUL     AX                ; q3*q3
  164.              NEG     AX                ; subtract
  165.              SBB     BX, DX            ;  q3*q3
  166.              SBB     SI, 0             ;   from remainder
  167.              SBB     BP, 0             ; BP:SI:BX:AX = remainder
  168.              POP     DX                ; DX:CX:DI = q1:q2:q3
  169.              JNS     $quot3_ok         ; remainder must be positive
  170. $corr_quot3: DEC     DI                ; try next smaller value for q3
  171.              STC                       ; adjust
  172.              ADC     AX, DI            ;  value
  173.              ADC     BX, CX            ;   of
  174.              ADC     SI, DX            ;    remainder
  175.              ADC     BP, 0             ;     according
  176.              ADD     AX, DI            ;      to
  177.              ADC     BX, CX            ;       changed
  178.              ADC     SI, DX            ;        value of
  179.              ADC     BP, 0             ;         quotient
  180.              JS      $corr_quot3       ; until remainder positive
  181. $quot3_ok:   XCHG    AX, DI            ; get result
  182.              MOV     BX, CX            ;  into DX:BX:AX
  183.              ADD     AX, 80h           ; round result
  184.              ADC     BX, 0             ;  (no tie case and
  185.              ADC     DX, 0             ;   no mantissa overflow possible)
  186.              POP     CX                ; get saved exponent
  187.              MOV     AL, CL            ; stuff exponent into result
  188.              AND     DH, 7Fh           ; clear sign bit
  189.              POP     BP                ; restore TURBO-Pascal frame pointer
  190.              RET                       ; done
  191. RSqrt        ENDP
  192.  
  193.              ALIGN   4
  194.  
  195. CODE         ENDS
  196.  
  197.              END
  198.